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I. INTRODUCTION 



The studies of entropic phase transitions in hard hypersphere systems in dimensions 
four, five, and higher, are at present an active field of research. The possibihty of precise 
quantitative studies appeared owing to the extension of molecular dynamics calculations to 
systems of hard hyperspheresi. Further development concerning the problem of freezing has 
been reviewed and discussed in^ where numerous references can be found. The motivation 
for studying fluids at D > 3 given in^ stresses the fact that the knowledge of fluid behavior in 
different dimensions can be effectively used to construct the corresponding density functional 
theory. Another reason is that in the theory of phase transitions one can expect important 
simplifications in dimensions D ^ 1. In the case of hyperspheres the solution ai D = oo 
could be used to develop a perturbative approach toward lower dimensions. 

The fluid to solid transitions in dimensions four, five, and six have been recently analyzed 
by advanced computations^,-,-. In paper^ both molecular dynamics and Monte Carlo sim- 
ulations have been used to study the onset of crystallization as reflected in the structure of 
the radial distribution function. Of course, the problem of primary importance is then the 
question of packing of hyperspheres. An interesting observation of geometrical frustration in 
four dimensions reported in^ is here a good illustration. The extension of numerical analysis 
of the fluid-solid transition to even higher dimensions six and seven can be found in^. 

The analysis of the fluid-crystal interfacial free energy in four, five , and six dimen- 
sions performed in^. showed that fluid stability increased with growing dimension. This 
interesting observation permitted to establish a connection with recent theories of jamming 
behaviour-,-. 

The instability of a hard hypersphere fluid with respect to a hypercubic crystal was 
analyzed in^. The authors mobilized density functional theory taking advantage of the 
exact relation for dimension D = 1 and for D = oo between the singlet density of an 
inhomogeneous system and the two-particle direct correlation function, and obtained via 
an analysis of bifurcations (seeii) an original estimate of the density of closest packing of 
hypercubic lattices. 

The question of phase structures appearing in hypersphere systems when the spatial di- 
mension D tends to infinity is the object of intensive studies. A general discussion of "magic 
dimensions" for which special lattice packings appear can be found However, it is 
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still by no means clear what kind of correlations persist when D ^ 1. The possibility of a 
simplification at D = oo has been strongly suggested by the study of Mayer series^ . The 
publications^,—, beyond reporting new results, provide a thorough description of the present 
state of the theory in high dimensions. A most interesting guess from existing results for- 
mulated ini^ is that in very high dimensions optimal packings of hard hyperspheres will be 
disordered, subject to decorrelation principle. This challenging hypothesis of disorder even- 
tually replacing closed packed crystalline structures became an important and fascinating 
theoretical question (see e.g. chapter 15.4 

Interpretation of growing wealth of precise numerical data requires a theory. Various 
theoretical approaches used for hyperspheres have been reviewed in^ where the density func- 
tional theory, virial expansions, scaled-particle theory, free-volume theory, Percus-Yevick 
and hypernetted chain integral equations are discussed. However, to our knowledge the 
superposition approximation, well known in the theory of liquids (for a critical review see^^) 
has not been systematically analyzed up to now in dimensions higher than three. The dis- 
cussion of the content of this theory for > 3 is our contribution to the current theoretical 
studies of hard hyperspheres. 

Our main object in the present paper is to investigate predictions of Kirkwood's superpo- 
sition approximation concerning the existence of phase transitions in hypersphere systems. 
The analytic and numerical results for D = 2, and D = 3, have been described in our 
previous work^^, showing good agreement with numerical studies. It turns out that we can 
apply the methods developed in=fi to investigate higher dimensions as well. In fact, we have 
at our disposal an analytic criterion for structural changes in arbitrary D. Comparison with 
existing numerical data confirms the correctness of predictions as far as the existence of 
structural transitions is concerned. So, although it is not known to what extent Kirkwood's 
approximation is valid for D ^ 1, the qualitative predictions in arbitrarily high dimensions 
are certainly worth examining. Owing to our method, we are also able to derive quantitative 
results, and compare them with numerical data. 

Our starting point is the equilibrium YBG hierarchy under Kirkwood's closure. In Section 
II the integral equation for the radial distribution function is derived. In Section III we solve 
the equation by iterations which permits one to determine the values of packing fractions 
corresponding to phase transitions in dimensions four, five, and six, and compare them with 
known numerical results. Section III is devoted to the discussion of — t- oo asymptotics. 
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The paper ends with concluding comments. 



II. SUPERPOSITION APPROXIMATION IN D DIMENSIONS 

The number density Us of s-particle configurations in which hard spheres of diameter 
a occupy space points (ri, r2, r^) can be conveniently written in terms of dimensionless 
positions = rj/a as 

n^(xi,X2, ...,x^) = n"" Jj6'(|xa - Xf,] - l)?/^(xi,X2, ...,xj (1) 

a<b 

where n is the number density of a uniform equilibrium state. The product of unit step 
functions 6'(|xa — X5I — 1) represents the excluded volume factor, and is the s-particle 
distribution, depending on dimensionless distances Xab = |xa — x^l, 1 < a < h < s. 

The second equilibrium Yvon-Born-Green (YBG) hierarchy equation for hard spheres in 
D > 2 dimensions has the form (see the derivation for D = 2 in^) 



-mix) = na^ j dcr{k ■ <t)6'(|x - <t| - l)yi,{x, 1, |x - cr\) (2) 



d 

dx'' 

Here x = |x|x = xx denotes the dimensionless relative position of a pair of hard spheres. 
Distances are measured in the sphere diameter a, so x = 1 describes a contact configuration. 
<T and X are unit vectors. The integration spreads over the solid angle with 

da- = sin'^~^(0i) sin^~'^(02)--- sm{(j)D^2)d(f)id(f)2---d(f)D-i (3) 
We choose the coordinate system such that 

X ■ <T = cos 01 

Under the Kirkwood superposition approximation the dimensionless three-particle den- 
sity ysi^x, 1, |x — <t|) factorizes into the product of two-particle distributions, and we get from 
02]) a closed nonlinear equation 
d 

dx' 

Introducing a simplified notation Y{x) = y2{x) we find that the correlation function H{x) 
Y{x) — 1 satisfies the equation 
d 



-mix) = na^ j d&{x ■ &)9{\x - &\ - I)y2{x)y2{l)y2i\x - ^\) (4) 



dx 



ln[H{x) + 1] = na^Y{l) j do-{x ■ o-)9{\x - 6-| - 1)[1 + H{\x - a] 



(5) 



We denote by v{l,D) the volume of a unit sphere in D dimensions 



T{l + D/2) 



(6) 



Using the formula 



2n 

2 I U(/^3--- / / d(f) D^iSm 

Jo Jo Jo 



we rewrite equation ([5]) in the form 



D-3i 



)...sm{<PD~2) = iD-l)v{l,D-l) (7) 



— \yi[H{x) + 1]=\{D) / dcj) cos sin^-2(0) e[x-2 cos 0] [1 + Hi^Jx"^ - 2x cos + 1)] (8) 



fix 
with 



A(D) =n(T^F(l)(D- 1) (9) 
The right hand side of equation (|8]) can be further simplified. It contains the integral 

e{2-x) 



cos sin (0)6'(x — 2 cos( 



D - 1 



1 - 



X 



21 {D-l)/2 



(10) 



and the integral involving the correlation function 

/•TT 

/ d(j) COS sin^^^(0)^^(x — 2 cos (f))H{-\/ x^ — 2x cos 0+1 
Jo 

"77/2 



(11) 



cos0sin^~^(0) e{x - 2 cos 0) H{^x'^ - 2xcos0 + 1 ) - H{\/x'^ + 2xcos0 + 1 ) 



Upon integrating both sides of equation ([5]) over the interval (x, oo) we get 



\n[H{x) + 1] = A(D) {Ri{x) + i?2(a;)} 



where 



i?i (x) 



D - 1 



rf2^(2 - z) 



z 

1-4 



(D-l)/2 



e{2 - x) 

D - 1 



1-^ 
4 



(D-l)/2 



(12) 



(13) 



and 



R2(X) 



/■oo r''^/'^ 

/ dz I (i0cos0sin 2(0)[/7(v/22 + 22COS0+ 1 )-^(^-2 cos 0) H{^z^-2z cos0+ 1 )] 
The term R2 can be transformed in the following way 

poo /'^/^ / 

R2{x) = J dz j rf0cos0sin^-2(0)[i/(^(;z + cos0)2 + sin2 0) (14) 



e{z -2 COS ^) H{^{z - COS + sin^ ) ] 

POO poo '\ I 

/ dz- dz 9 {z - cos (f)) [ H{^z^ + sm\(j))) 

Jx+costp Jx—cos<f> ) 

/ d0cos0sin^~2(0) / dz 9 {z - cos (fy H{^J z'^ + sin^icj))) 

Jo Jx—cos<f> 



Tv/2 

d(f) COS (f) sin^~^(0) 



Introducing the integration variable /i = sin we get 

p1 r+oo 

R^[x) = - d^iii^-^ dz9{z'^ + - l)9{^l - - \z - x\) z"^ + y?) (15) 
Putting now s = \/ z^ ^ leads to the equality 

Mx) ^-Jds sH(s)e{s -1)1^1' J^'l/ ''^^^'- ~ ' ' 

^ - / - 1, [■ - £l±£l^) 

One can further simplify this expression by using the integration variable w = \/\ — v^. 
Indeed, we have 

\^'' du -4^9 ( - ^!±f!^^ = dw (l+^2)(D-3)/2^ f _ X' + S'-1 \ 

Jo V 2sa; J J^y, ^ ' \ 2xs J 

(17) 

= r rf^^^(l_|a;_s|)(l_^2)(D-3)/2 
J {x^+s2-l)/2xs 

The last equality follows from the fact that 



(x2 + S^-1) yi^^ 
- > 



2xs 

Using this result we eventually find 

px+l nl 

R2ix) = - dss''-'H{s)9is-l) rfw(l-M;2)(^-3)/2 (18) 

Jx-l J {x^+s'^-l)/2xs 

We can thus write the integral equation ( IT2l) for the two-particle correlation function H{x) 
of D— dimensional hypersheres in the form 

H{x) = CH{x) = -1 + exp{X{D)[Ri{x) + R2{x)]} (19) 

where 



Ri{x) + R2ix) = { ^"^^ I dwil-w'Y''-'^/^ 

lx/2 



f 29{2-x) f 
\ D-1 7, 

dss''-^H{s)9{s-l) f rfw(l-«;2)(D-3)/2 

Jx-l J (x^+s^~l)/2xs 



III. ITERATIVE SOLUTION OF THE INTEGRAL EQUATION 



The integral equation (fT9|) is solved by a standard Neumann method with succesive over- 
relaxation. The iterative solutions are then given by 

Hn = il- a)Hn-i + aC{Hn-i) (20) 

where C has been defined in equation f|T9|) . The relaxation parameter a was taken to be 0.1 
(for D = 3, 4), and 0.05 (for D = 5, 6). Iterations were continued until successive values 
of H{x = 0) differed by less than e = 10~^, except in the vicinity of the threshold volume 
fraction 0* (see below), where the convergence was slow and the iterations were discontinued 
at e = 10-2. 

Examples of correlation functions obtained in this way are presented in Figs.d], [2J Clearly, 
with increasing volume fraction the decay of H{x) becomes slower, and a pronounced peak 
structure appears. 

The comparison of our results with the molecular dynamics simulation data of Estrada 
and Robles^ presented in Fig.Oshows that the correlation function obtained from the integral 
equation fll9p has a lower contact value, H{1), and shifted maxima with respect to the 
molecular dynamics curves. Analogous quantitative differences between the predictions of 
Kirkwood's approximation and the molecular dynamics data are also observed at higher 
dimensions and other volume fractions. 

In order to investigate the possibility of a structural change we proceed as in^*^ by con- 
sidering the form of equation for large distances x. Using the asymptotic formulae 
ln.[H{x) + 1] = H{x), and \/ x^ — 2xcos{iIj) + 1 = x — cos{iIj), we get an integral equation 

^H{x) = \{D)[ #cos^ sin^-2(V^)iJ(x-cosV^) (21) 
dx J Q 

We then consider the solution H{x) as a linear combination of exponential modes 

H^,{x) = exp{K,x), 
where k is a complex number satisfying the equation 

/•TT 

K = X{D) / dijj cos ijj sm^~'^ (ijj) exp[—K cos ijj] (22) 
Jo 

We look for the mode with the slowest decay. The disappearance of exponential damping in 
this mode announces the change in the nature of correlations, and thus implies a structural 
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FIG. 1. Pair correlation function, H{x), for the hard-sphere gas {D = 3) at the volume fraction 
(p = 0.40 (dashed line) and (f) = 0.51 (solid line). 

change. Such a possibility is equivalent to the appearance of a purely imaginary solution 
K = ib of fl22|) . with b obeying 



D/2 



JD/2{b) = 



In deriving equation (!23|) from (!22|) we used the relation 

'2 



n \ - j T{iy + 1/2) J„{z) = J dijj sin^" {^p) cos{z cos ijj) 



together with 



dz z"" 



Ju+l 



Equation (123|) has a solution if and only if X{D) > X*{D) where 

■JD/2ib) 



(23) 



(24) 



(25) 



(26) 



In order to evaluate the absolute minimum in fl26|) we note that according to relation f l25|) 
all extrema of function Jy(b)/b'^ for 6 7^ are attained at points which are zeros of function 
Ju+iip). One finds 



Min 



JD/2ib) 



Jd/2[j(1 + ^/2,1)] 



(27) 



[j(l + /^/2,l)p/2 

where j(l + D/2, 1) is the first positive zero of function Ji+d/2- Inserting (1271) into (126!) we 
find the following values of \*{D) : 
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FIG. 2. Pair correlation function, H{x), for the hard-hypersphere gas {D = 5) for the volume 
fraction (J) = 0.10 (dashed line) and (j) = 0.27 (solid line) . 




FIG. 3. Pair correlation function, H{x), for the hard-sphere gas {D = 3) at the volume fraction 
(j) = 0.40 obtained from the solution of the integral equation (jl2p (solid line) compared with the 
molecular dynamics simulations of Estrada and Robles^ (dashed line) . 



9 



109 
10^ 
10^ 







20 



40 



60 



80 



D 
100 



FIG. 4. Critical value A* as function of D. 



A*(3) = 17.407 
A* (4) = 43.44 
A*(5) = 91.23 
A* (6) = 172.76 

We note here a rapid increase of the critical value X*{D) with dimension. Indeed, the formula 
(I26p implies a rapid growth (faster than (e/2)-^/^) illustrated on Fig. HI (In our paper— we 
found at D = 3 the threshold value 34.81 for the parameter considered by Kirkwood. et 
al.— , and equal to 2A of the present paper). 

The function \{D) defined in (j9]) is related to the volume fraction 4>{D) occupied by the 
spheres 

<l>iD) = n{£l''vil,D) (28) 



by the equation 



where 



X{D) = 2^y(l)(D - 1)%^^^0(D) 



v(l,D-l) 



v{l,D) 

r(i + D/2) 



v{l,D) v^r[(l + D)/2] 
10 



(29) 



(30) 



Of course, the contact value of the radial distribution also depends on the volume fraction 
which we will highlight in what follows by the notation Y{l,(f)). Example of the resulting 
A(0) dependence (for D = 4) is given in Fig. O 

Taking from the iteration procedure we can estimate the hyper-volume fractions 

4>*{D) corresponding to the above-derived critical X*{D). In this way, we get 

0*(3) = 0.52 
0*(4) = 0.40 
0*(5) = 0.28 
0*(6) = 0.21 

The contact values Y{l,(j){D)) evaluated at 4'*{D) show with increasing D a decreasing 
behavior 

F(l,0*(3)) = 2.79 

F(l,0*(4)) = 2.67 

y(l,0*(5)) = 2.61 

F(l,0*(6)) = 2.51 

The Kirkwood approximation does predict phase transitions in dimensions 3, 4, 5, and 6 
because the threshold values of the volume fractions given above are lower than the optimal 
volume fractions (f)max{D). Indeed, in three dimensions we know the exact result 0maa(3) = 
71 / 3^/2 = 0.7404, whereas in dimensions 4 ,5, and 6 the largely accepted conjectures for the 
densest lattice packings^,^^ yield 0™ax(4) = 7rVl6 = 0.6168, 0™ax(5) = 27rV30v/2 = 0.4652, 
0™ax(6) = 7rV48y3 = 0.372. 

The numerical results reported in^ show that the phase coexistence region for D = 3 is 
in the density range 0.494 < < 0.54, whereas for D = 4 and D = 5 the corresponding 
ranges are given by 0.288 < (p < 0.337 and 0.174 < (p < 0.206, respectively. Finally, 
for D = 6 the estimated coexistence range corresponds to 0.105 < < 0.138. Thus the 
Kirkwood approximation overestimates the transition point, situating it for D = 4,5, and 
6 beyond the coexistence region. The discrepancy between our results and the numerical 
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ones measured by the ratio \(p*{D) — (f)fr{D)\/(j)fr{D), where (f)fr{D) is the volume fraction 
at freezing, increases with growing dimension. Moreover, the numerical simulations of Ref.- 
place the volume fraction of maximally random jammed state at = 0.2 ± 0.01 for D = 6, 
which suggests that the transition at (f)*{D = 6) = 0.21 might be kinetically inaccessible. 

As far as the origin of the above discrepancies is concerned the following remark can 
be made. According to a thorough and subtle analysis of three-particle correlations in 
hard spheres performed by B. J. Alder— the superposition approximation gives very good 
quantitative results provided one extracts the radial distribution directly from the triplet 
distribution without using the YBG hierarchy. His important conclusion is that poor quan- 
titative results of the Kirkwood approximation "...are due to an extreme magnification of 
the error by the integral equations in which it was introduced." So, it seems possible that 
the above mentioned discrepancies between our results and numerical predictions are mainly 
due to the fact that we apply the superposition approximation to the hierarchy equation. 
However, the qualitative conclusions may be correct. 

At this point another important problem must also be considered. The discussion of 
our results presented so far assumes that the threshold volume fraction (f)*{D) corresponds 
to crystallization. But this cannot be really inferred from our approach. All we know 
is that for (f){D) > 4'*{D) correlations change their nature, and the law of exponentially 
damped oscillations must be replaced by another one. The new law must describe long-range, 
nonintegrable correlations. This last condition is obviously satisfied by crystal structures. 
But it would be also satisfied by states with correlations decaying according to nonintegrable 
power laws. To our knowledge, such power laws have not been found in jammed or glassy 
states^i^. Consequently, the Kirkwood threshold volume fractions 4'*{D) cannot correspond 
to their appearance. We shall come back to this question at the end of the next section. 

IV. D^oo: POSSIBILITY OF STRUCTURAL CHANGES 

We now turn to the discussion of the high-dimensional asymptotics Z) — )■ cxd having 
in view the answer to the fundamental question of existence of phase transitions within 
Kirkwood's approximation. It follows from the analysis presented so far that correlations 
change their nature provided X{D) > X*{D). We have thus to investigate the content of the 
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FIG. 5. The dependence \{(t)) for D = 4. The dashed hne indicates the threshold value, A* = 43.44. 



inequality 



X{D) > \*{D) 



21-D/2 



[j(l + /^/2,l)]^/2 



v/^r[(D - l)/2] {Jn/2{3{l + D/2,l) 
Using the relation 

Mj{l + u,l)] = Jl^,[jiu + l,l)] 
together with asymptotic formulae for large u (see2^) 

3i/6r(2/3) 



, 1) ^ zy + const z/^/^ Jli 



we find 



I3>1 



2V3vrz/2/3 
21-D/2 (1 + D/2)^/22i/3 



V^T[{D - l)/2] 3i/6r(2/3) 
In view of relations fl29l) and fl30l) we can rewrite fl34l) as 



1 + /^/2) 



2/3 



2^{D - l)y(l,0*(Z})) i^//"'"^(f . 0*(Z}) 



r((i + z})/2) 

21-D/2 (1 + £)/2)^/22V3 



1 + /^/2) 



2/3 



(31) 
(32) 

(33) 
(34) 

(35) 



T[{D - l)/2] 3i/6r(2/3) 
A straightforward calculation yields then the following large D formula for the threshold 

volume fraction (j)*{D) 



cf>*{D)Y{l,<p*{D))= const (f ) (^) 



e 



(36) 
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where 

const = 2^/3g/[y2^3i/6p(^2/3) 

In order to check whether the volume fraction (f)*{D) satisfying (136|) can be attained one 
needs a precise knowledge of the upper bound for possible volume fractions in D dimensions. 
One also needs the behavior of the contact value Y{1, 0(-D)) for D — t- oo. 

According to the Mayer series study^^ the hard-hypersphere equation of state at D = oo 
has a remarkably simple form 

p = nkBT[l + ^na^v{l, D)] (37) 

whereas the exact equation reads 

p = nkBT[l + ^na^'vil, D)Y{1, 0(D))] (38) 

It follows that Y{1, 4>{D)) = 1 for D = oo. In fact, we have noticed within our approach the 
decrease of Y{l,(j){D)) with growing dimension when analyzing data for D = 3,4,5 and 6. 
This observation suggests that Y{l,(f){D)) could monotonously approach 1 when D — > oo. 
We will thus consider Y{l,(f)*{D)) in ( l36l) for D ^ 1 as a number close to 1. The resulting 
scaling of the volume fraction at phase transition reads 



1/6 . g D/2 



(39) 



A straightforward calculation shows that the Rogers rigorous upper bound for lattice 
packings^,— 

^(D) < ^ (40) 
is satisfied by scaling ( 139|) of (p*{D). Indeed 

2^/2e 
lim 0*(D) = 

Also the stronger Kabatiansky and Levenshtein bound^^ 



20.59901) 



does not lead to contradiction with ( 139|) for D — )■ oo. We thus conclude that scaling ( 139|) is 
compatible with existing upper bounds for crystals. 
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The crossing of the threshold volume fraction (p*{D) leads to a change in the structure 
of hyperspheres reflected by the appearance of long range correlations. As we have already 
remarked at the end of Section 3, this is the reason why we do not expect (f)*{D) to announce 
the passage to glassy or jammed states. Let us just note that the scaling 

(f)K{D) = 2-^ DlnD (41) 

of the Kauzmann point (pxiD) of the thermodynamic glass transition^,— follows quite a 
different law compared to that derived for 0*(-D) in (!39|) . We note that 
lim£)_!.oo 0_ft'(-D)/0*(D) = which implies that for sufficiently high D the glass transition 
would occur at lower volume fraction than the transition at (f)*{D). 

V. CONCLUDING COMMENTS 

Our main object in this paper was the investigation of the possibility of structural tran- 
sitions in hypersphere systems within Kirkwood's superposition approximation. To this end 
we employed the simple criterion derived from the equilibrium YBG hierarchy: exponential 
damping of the oscillating pair correlation function H{x) disappears when the dimensionless 
parameter X{D) = ncr^{D — 1)Y{1, (f){D))v{l, D — 1) attains the threshold value X*{D) given 
by equation (d^. For X{D) > X*{D), the large distance behavior of correlations is necessar- 
ily changed. In order to check whether the transition is possible we had to make sure that 
the threshold volume fraction (f)*{D) corresponding to X*{D) was smaller than the maximal 
possible value (pmax- This was the most difficult point because it required the knowledge of 
the contact value of the radial distribution Y{l,(f)*{D)), and thus the solution of the inte- 
gral equation (fT9|) . We performed this program for D = 3,4,5, and 6 concluding that the 
superposition approximation does predict phase transitions for 3 < < 6, in accordance 
with numerical results, and is thus qualitatively correct. However, it yields threshold values 
of the volume fraction higher than those following from numerical studies for crystallization. 

Our investigation of the situation at D = oo permitted to derive the asymptotic form of 
scaling fl5I?l) for the volume fraction at phase transition showing consistency of the superpo- 
sition approximation with crystallization. We checked that the known lattice upper bounds 
for the maximal volume fractions in D dimensions are not restrictive enough to eliminate 
the possibility of crystallization in arbitrarily high dimension. 
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An interesting question left open is the limit lini/j^oo ^(1) 4'*{D)). We noticed in Section 
III that the contact values Y{l,(j)*{D)) decreased with increasing dimension D. According 
tcM; the contact value at = oo is simply equal to 1. The evaluation of the above limit 
within the superposition approximation would be thus an important test for this theory. 

Let us finally stress the fact that although our approach predicts disappearance of the 
fluid structure characterized by exponentially damped correlations, it cannot predict the 
precise nature of the new emerging phase. Clearly further work needs to be done on this 
point. 
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